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Abstract. Slow-fast dynamical systems have two time scales and an explicit parameter repre- 
senting the ratio of these time scales. Locally invariant slow manifolds along which motion occurs 
on the slow time scale are a prominent feature of slow-fast systems. This paper introduces a rig- 
orous numerical method to compute enclosures of the slow manifold of a slow-fast system with 
one fast and two slow variables. A triangulated first order approximation to the two dimensional 
jy^ | invariant manifold is computed "algebraically". Two translations of the computed manifold in the 

fast direction that are transverse to the vector field are computed as the boundaries of an initial 
enclosure. The enclosures are refined to bring them closer to each other by moving vertices of 
the enclosure boundaries one at a time. As an application we use it to prove the existence of 
tangencies of invariant manifolds in the problem of singular Hopf bifurcation and to give bounds 
on the location of one such tangency. 
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1. Introduction 

Invariant manifolds and their intersections are important features that organize qualitative 
properties of dynamical systems. Three types of manifolds have been prominent in the subject: (1) 
compact invariant tori [23], (2) stable and unstable manifolds of equilibria and periodic orbits [2~T| 
[7], and (3) slow manifolds of multiple time scale systems [18J. Interval arithmetic and verified 
computing have been used extensively to give rigorous estimates and existence proofs for invariant 
tori and occasionally to locate stable and unstable manifolds, but this paper is the first to employ 
these methods to locate slow manifolds. Each of these three cases pose numerical challenges to 
locate the manifolds. 
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Many methods that locate invariant tori assume that the flow on the tori is smoothly conjugate 
to a constant flow with dense orbits. Existence of this conjugacy confronts well known small 
divisor problems and the winding vector of the flow must satisfy diophantine conditions in order 
for this problem to be solvable. Typically, the numerical methods produce a Fourier expansion of 
the conjugacy which is determined up to a translation. The manifolds are located by projection 
onto a discrete set of Fourier modes and solving a fixed point equation for the coefficients of the 
conjugacy. 

The computation of stable and unstable manifolds of equilibria and periodic orbits is a "one- 
sided" boundary value problem. The manifolds consist of trajectories that are asymptotic to the 
equilibrium or periodic orbit. In the case of an equilibrium point of an analytic vector field, the lo- 
cal stable and unstable manifolds are analytic graphs that have convergent asymptotic expansions 
whose coefficients can be determined iteratively. The most challenging aspect of computations 
of two dimensional manifolds arises from the way that trajectories do or do not spread out in 
the manifold as one departs from the equilibrium or periodic orbit. As illustrated by the Lorenz 
manifold [21J, the manifolds can twist and fold in ways that present additional geometrical compli- 
cations for numerical methods. The development of rigorous bounds for these invariant manifolds 
follows similar principles to the verified computation of individual trajectories. 

Multiple time scale vector fields, also known as singularly perturbed differential equations, occur 
in many settings: systems of chemical reactions, lasers, fluid dynamics and models of the electrical 
activity of neurons are a few examples. Borrowing terminology from fluid dynamics, the solutions 
of these systems can have (boundary) layers in which the fast time scale determines the rate at 
which the solution varies as well as long periods of time during which the solution evolves on the 
slow time scale. The slow motion typically occurs along slow manifolds that are locally invariant. 
The slow manifolds play a prominent role in qualitative analysis of the dynamics and bifurcations 
of multiple time scale systems. Indeed, model reduction procedures are frequently employed that 
replace a model by a lower dimensional model that approximates the motion along a slow manifold 
and ignores the fast dynamics of the original model. The ideal for this type of model reduction is 
an algorithm that computes the slow manifold exactly. That ideal seems very difficult to achieve 
and is not addressed in this paper. Instead, we seek rigorous bounds for the location of the slow 
manifold that are tight enough to give information that can be used in the analysis of bifurcations 
of the system. 

To explain the methods we introduce in the simplest terms, we focus upon slow-fast systems 
that contain an explicit parameter e that represents the ratio of time scales. Moreover, we restrict 
attention to systems that have two slow variables and one fast variable and use a single example as 
a test case. In principle, the methods generalize to the case of codimension one slow manifolds, and 
the definitions and existence proofs in Sections [2] and [3] have obvious higher dimensional analogues. 
In practice, however, due to the scarcity of tools for computational geometry in higher dimensions, 
implementing a higher dimensional version would be a significant extension of the work described 
in this paper. We comment on generalizations from the setting of systems with two slow and one 
fast variable in the discussion at the end of the paper, but leave consideration of further details 
to future work. 

Slow manifolds of multiple time scale systems present unique theoretical and numerical chal- 
lenges compared to the computation of invariant tori and (un) stable manifolds. The first of these 
challenges is that theory is developed primarily in terms of "small enough" values of the parameter 
e measuring the time scale ratio of a slow-fast system. Numerically, one always works with specific 
values of e. The convergence of trajectories as e — > is singular, making it difficult to develop 
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methods framed in terms of asymptotic expansions in e. Divergent series are the rule rather than 
the exception in this context. The rich history of numerical integration methods for stiff systems 
and the large literature on reduction methods for kinetic equations of chemical systems reflect the 
difficulty of computing attracting slow manifolds, the simplest case for this problem. Comput- 
ing slow manifolds of saddle-type presents the additional challenge that most nearby trajectories 
diverge from the slow manifold on the fast time scale in both forward and backward time. The 
second theoretical difficulty in finding slow manifolds is that they are only locally invariant in 
most problems of interest. The local invariance is accompanied by a lack of uniqueness: possible 
manifolds intersect fast subspaces in open sets whose diameter is exponentially small in e; i.e., 
bounded by exp(— c/e) for a suitable c > 0. Methods based upon root finding of a discretized set 
of equations must choose a specific solution of the discretized equations. 

We compute enclosures of slow manifolds by exploiting transversality properties that improve 
as e — > while being suitable for fixed values of e. The methods do not identify a unique object 
and are well suited to locating locally invariant slow manifolds. If H is a hypersurface and F is 
a vector field, then transversality of F to H is a local property: verification does not rely upon 
computation of trajectories of F. For a slow-fast vector field with one fast variable, translation of 
a normally hyperbolic critical manifold along the fast direction produces a transverse hypersurface 
when the translation distance is large enough. Translation distances proportional to e suffice. In 
this paper, we use piecewise linear surfaces H as enclosing manifolds. For the example we consider, 
transversality at vertices of a face of H implies transversality of the entire face. This reduces the 
computational complexity of checking transversality sufficiently that iterative refinement of the 
enclosures was feasible. 

Since slow manifolds are objects that are defined asymptotically in terms of the parameter e, 
they are not directly computable using finite information. One part of this paper is devoted to 
the development of a mathematical framework within which slow manifolds are defined for fixed 
values of e > 0. We define computable slow manifolds and relate this concept to the slow manifolds 
studied in geometric singular perturbation theory. All computations and statements in this paper 
are for computable slow manifolds. This is similar in spirit to the finite resolution dynamics 
approach of Luzzatto and Pilarczyk |26) . 

Our work is motivated by the study of tangencies of invariant manifolds. Significant global 
changes in the dynamics of a system have been observed to occur at bifurcations involving tan- 
gencies. Proving the existence of tangencies is intrinsically complicated because the manifolds 
themselves must be tracked over a range of parameters. Computer-aided proofs of tangencies of 
invariant manifolds have previously been studied by Arai and Mischaikow in [Tj, and Wilczak and 
Zgliczyfiski in [38] . In Section [TJ we prove that a tangency bifurcation involving a computable 
slow manifold occurs in the singular Hopf normal form introduced in [§]. 

1.1. Slow- fast systems. Slow-fast differential equations have the form: 

(1) ex = f(x,y,e) 

y = g(x,y,e), 

where i 6 R" j 6 R"\ f : R"+m+i _^ R » and g . R n+m+i _^ R m_ We assume t h at t h e vector 
field (f,g) is smooth (C°°), although most of this paper can easily be adapted to the finitely 
differentiable setting. Here x and y are the fast and slow variables, respectively. Throughout the 
paper we consider the case m = 2 and n — 1 of two slow variables and one fast variable. 
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We define the critical manifold, as the set 

(2) S :={(x,y)eR n+m :f(x,y,0) = 0}. 

The critical manifold is normally hyperbolic at points where D x f is hyperbolic; i.e., has no eigen- 
value whose real part is zero. Points where So is singular are referred to as folds. On the normally 
hyperbolic pieces of the critical manifold, x is given as a function of y, x = ho(y). The correspond- 
ing differential equation 

(3) y = g(h o (y),y,0) 

is called the slow flow. If one instead rescales time with e and puts e — in (JT]), one gets the layer 
equation: 

(4) x' = f(x,y,0) 

V = 0, 

Note that the manifold So is exactly the set of critical points for the layer equation. Singular 
perturbation theory studies how the solutions to ([T]) for e small, but positive, can be understood 
by studying solutions to ([3]) and (j4j). 

When S*o is normally hyperbolic and e > is sufficiently small, geometric singular perturbation 
theory [18] ensures that the critical manifold perturbs to a slow manifold. Slow manifolds are 
locally invariant and 0(e) close to the critical manifold. However, slow manifolds are not unique, 
although different choices are within 0(e~ c ^ £ ) distance from each other. We denote slow manifolds 
by S^.The purpose of this work is to compute approximations of S E that are guaranteed to be 
of a certain accuracy. This is achieved by computing two approximations that enclose the slow 
manifold. The two approximations of the slow manifold are triangulated surfaces transverse to the 
vector field. To prove the transversality, we use interval analysis, to be explained in Subsection [T72J 
Interval analysis is a general technique that enables mathematically rigorous proofs of inequalities 
on a digital computer. 

To simplify notation we denote the two slow variables by y and z, i.e., from now on y £ R, and 
the vector field in the slow variables is denoted by g = (g y , g z ). We also assume that /, g y , and g z 
are independent of e. To summarize, the systems we study are of the following form: 

(5) ex = f(x,y,z) 

V = 9 v {x,y,z) 
z = g z (x,y,z), 

where x,y, z £ R, and /, g y , g z : R 3 — > R. We will sometimes use the notation F = (f, g y , g z ). 

1.2. Validated numerics. Interval analysis was introduced by Moore in [27] as a method to use 
a digital computer to produce mathematically rigorous results with only approximate arithmetic. 
Tucker [3l] is a modern introduction to the subject, and more advanced topics are discussed by 
Neumaier [31] . The main idea is to replace floating point arithmetic with a set arithmetic; the 
basic objects are intervals of points rather than individual points. Together with directed rounding 
this method yields an enclosure arithmetic that allows for the rigorous verification of inequalities. 
To use interval analysis to produce a mathematical proof, often called ( auto- jvalidated numerical 
methods, one has to prove that the statement at hand can be reduced to a finite number of 
inequalities, and then verify that these inequalities are satisfied. Interval arithmetic is used for the 
verification. The objects used to describe sets in validated numerics are typically convex sets in 
some coordinate system, e.g., intervals, parallelograms, or ellipsoids. In this paper we will employ 
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triangular meshes of surfaces, an approach that previously, in this setting, only has been used in 

1.3. Computation of invariant manifolds. The study of invariant manifolds [12] is central 
to the theory of dynamical systems. The behavior of a system can often be understood by un- 
derstanding its invariant structures. Numerical computations of invariant manifolds El IH El 
[TO] [T5| l21j |32| |33| 136] are important in many applications. There are no universally applicable 
methods to compute invariant manifolds; to be efficient, they have to be tailored for the specific 
class of problems one is studying. Computing invariant manifolds of slow-fast systems is partic- 
ularly challenging. Two existing methods are [10], and no rigorous methods exist. The main 
idea of our method is to refine a first order approximation of the manifold by local modifications 
that maintain transversality of the enclosing manifolds. Interval arithmetic is used to make the 
local computation of transversality rigorous. This is similar in spirit to the methods developed 
in [5] to study the phase portraits of planar polynomial vector fields. Even in the planar case 
the verified computation of phase portraits is a challenging task, and the few methods that exist 
include [H HH]. 

2. Overview of the method 

This section describes our method to compute enclosures of the slow manifold of a slow-fast 
system of the form (J5J). We start by giving an overview of the main ideas of the method. There 
are five main steps in the algorithm: 

(1) triangulation of the critical manifold, 

(2) computing the 0(e) correction term for the slow manifold, 

(3) constructing left and right perturbations of the slow manifold, 

(4) proving that the left and right perturbations enclose the manifold, and 

(5) tightening the enclosure by contracting the left and right perturbations towards each other. 

The first step is to compute a triangulation of the critical manifold, which is adapted to its 
geometry. The manifold is defined implicitly by the condition f(x,y,z) = 0. In the example 
we consider in Section [51 we solve this equation to obtain explicit expressions for the functions 
of the form x = ho{y,z) whose graphs lie in the critical manifold. Alternatively, one computes 
approximations to ho using, e.g., automatic differentiation and continuation procedures. There are 
many software packages to compute triangulations of surfaces; we use CGAL |30] via its matlab 
interface. When a part of the critical manifold is represented as the graph of a function h , its 
domain in the plane of the slow variables can be triangulated, and then this triangulation can be 
lifted to the graph, as illustrated in Figure [TJ So that the triangles in the lifted triangulation have 
similar diameters, we choose triangles in the plane of the slow variables to have diameters that 
depend upon the gradient of ho- We stress that the rest of the algorithm is independent from how 
the triangulation of the critical manifold is constructed. Rather than using axis parallel patches, 
one could, e.g., use approximate trajectory segments of the reduced system to determine the piece 
of the domain of the slow variables, where the slow manifold is computed. 

We compute an approximation to the slow manifold using a procedure similar to that employed 
in stiff integrators that use Rosenbrock methods [13] . The tangent space to the critical manifold 
is orthogonal to the vector df. According to the Fenichel theory, the slow manifold is 0(e) close to 
the critical manifold in the C l topology, so its tangent space is approximately normal to df. At a 
point (x, y, z) in the (lifted) triangulation of So, we look for a nearby point (x + 5, y, z) at which the 
vector field is orthogonal to df(x,y,z). Since f(x,y,z) = and the normal hyperbolicity implies 
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(a) y (b) 

FIGURE 1 . The mesh generated for the example in Section [5j There is a fold at 
{y = 0}. (a) The Delaunay triangulation of the (y, z) plane that is generated by the 
geometry adapted mesh points, (b) The lift of the triangulation to the critical 
manifold. 

that dj ± 0, 

r (d y fg y + d z fg z ) 
£ (dj) 2 

is an approximate solution to this equation. Setting 5 to this value, we take (x+5, y, z) as a point of 
the triangulation of the approximate slow manifold. The critical manifold and the approximation 
to the slow manifold are illustrated in Figure [2] (a) and (b), respectively. We next perturb this 
triangulation of the approximate slow manifold in both directions parallel to the x-axis, as in 
Figure [2] (c), by a factor 2 J_6 5, where j is a natural number that will be specified later. In case 
that 5 is very small, we replace it by a 0(e 2 ) term. This procedure yields two surfaces that are 
candidates for the enclosing surfaces that we seek. 




FIGURE 2. Construction of the enclosing triangulations. This figure shows the 
projection on (x, y) coordinates of: (a) the critical manifold (solid) (b) the critical 
manifold (solid), and the slow manifold (dotted) (c) the critical manifold (solid), the 
slow manifold (dotted), and the two enclosing surfaces (dashed). 



To verify that the surfaces enclose the slow manifold, we check whether the flow of the full 
system ([5]) is transversal to the candidate surfaces. As the candidate surfaces are piecewise linear, 
we have to define what we mean by transversality at the edges and vertices of the triangulation. 
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Definition 2.1. Let Tel 3 be a triangulated, piecewise linear two dimensional manifold T = 
IjTj. Since T is a manifold, it locally separates R 3 into two sides. We say that a vector v is 
transverse to T if v and —v point to opposite sides ofT. A smooth vector field is transverse to T 
if it is transverse to T at every point of T . 

Figure [3^ a) illustrates this definition. Trajectories of the flow generated by F will all cross T 
from one side to another if F is transverse to T. If 71 and J~2 are triangulated surfaces transverse to 
the flow with opposite crossing directions, then they form enclosing surfaces for the slow manifold 
we seek. 



T 




FIGURE 3. (a) Transversality check on an edge of the triangulation. F is transversal 
if it does not belong to the cone C. (b) Transversality check on one face of the 
triangulation. To verify that the flow intersects the surface transversally, it suffices 
to prove that the vector field is never orthogonal to the normal of the surface, which 
a constant vector. So, we compute the range of F ■ n, and prove that it does not 
contain 0. 

Transversality is a condition that is local to each face of the triangulation, so we can check it 
on each face of the triangulation separately. To check the transversality condition on one face, 
we estimate the range of the inner product of the vector field with the normal of the face, as 
illustrated in Figure [3](b). Details about the existence of locally invariant, normally hyperbolic 
manifolds inside the enclosure are addressed in Section [3] below. 

The final part of the algorithm is to iteratively update the location of the vertices by moving 
them towards each other in small steps along the fast direction. We check that the transversality 
properties still hold, see Figure HI This tightening step is stopped when no more vertices can be 
moved. Note that the vertices of all 4 triangulations: the critical manifold, the approximate slow 
manifold, and the two perturbed manifolds, all have the same (y, z) components. 

3. Existence of locally invariant manifolds 

The method outlined in the previous section constructs two triangulated surfaces, in the phase 
space of a slow-fast system, that are transversal to the flow for the given e. In this section we 
discuss the existence of locally invariant manifolds enclosed between these two triangulations. We 
denote the two enclosing surfaces by £ and TZ, and the region enclosed between them by C. Note 
that L and TZ are graphs over the same compact region, so C is well defined. Specifically, if for 
some compact set of slow variables D, C = {(x,y,z) : x = hi(y,z), (y,z) G D}, TZ = {(x,y,z) : 
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FIGURE 4. Updating the enclosures of the slow manifold. This figure show the 
projection on (x, y) coordinates of: (a) The initial enclosures, (b) The enclosures 
are updated vertex wise, here the first half of the vertices are updated, (c) The new 
enclosure. 



x = h r (y,z),(y,z) e D}, then C = {(x,y,z) : x e [hi(y,z),h r (y,z)],(y,z) e D}. We would 
like to claim that there is a locally invariant manifold inside of C which is a graph over the slow 
variables. We must thus verify that it is possible to choose a subset of C which is a C l manifold, 
locally invariant, and whose projection onto the domain in the slow variables is bijective. We 
start this section by defining computable slow manifolds as objects associated to a fixed e. Similar 
to a slow-manifold, a computable slow manifold, is not unique. Informally, a computable slow 
manifold is a manifold close to the critical manifold where the flow is slow. We measure slowness 
by comparing the slopes of trajectories within our enclosure with the slope of the critical manifold. 
The relative slope is defined as a bound on the slope of trajectories divided by the slope of the 
critical manifold: 
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Definition 3.1. A computable slow manifold is a C 1 locally invariant, normally hyperbolic man- 
ifold, of the same dimension as the critical manifold, which projects injectively along the fast 
variable to the critical manifold, such that the its relative slope satisfies s(e) < 4|. 

Note that the order of s(e) is O away from the slow manifold, so the definition is consistent 
with the standard perturbative definition of slow manifolds [I8j . Slow manifolds are widely used 
in studies of slow-fast systems arising from biological or chemical models. However, computable 
slow manifolds are often the objects that are identified in applications: a locally invariant manifold 
at a fixed value of epsilon that follows the critical manifold closely, and on which the flow is slow 
[6j [10J HI]. This concept is captured by the definition of the computable slow manifold. Thus, 
our enclosures method gives a general and robust method to compute where candidates for such 
manifolds might lie in the phase space. 

We will explain why computable slow manifolds exist within C in the following special case, 
which is sufficient for the purpose of this paper. 

Assumption 3.2. Assume that: 

i All trajectories of C reach its boundary in forward and backward time. 
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ii The boundary of C, dC is piecewise smooth. Tangencies of the vector field with dC are 
quadratic (i.e., folds in the sense of singularity theory), and these tangencies occur along 
smooth curves that connect £ and TZ. 

iii There are invariant horizontal and vertical cone fields on C, and the vertical invariant cone 
field contains the fast direction of the vector field on C. 

Assume that the vector field is inward on C and TZ and denote by C; n and C out the sets in 
dC — C — TZ where the vector field points inward and outward, respectively. Choose a smooth 
curve, x = r (y,z), in dC — C — TZ such that the projection of the curve to the slow variables 
contains the projection of Ci n to the slow variables, and points on the curve on C ou t are images 
of the flow of points on the curve on Ci n . Flow this graph forward until each trajectory leaves C. 
The set swept out by these trajectory segments is: 

(7) S e := {0 t (x, y,z) : x = r (y, z), t (x, y, z) G C}. 

The set S e is well-defined, as smooth as ro and the vector field, and diffeomorphic to its projection 
onto the critical manifold So- Inflowing trajectories of C must exit through C ou t- The exit time is 
uniformly bounded, since C is compact. Hence, S e is well defined. The existence of invariant cone 
fields, with a normal vertical cone field containing the fast direction, ensures that S £ is a graph 
over the slow variables, and thus diffeomorphic to the corresponding part of the critical manifold. 
The final requirement of the definition - that the relative slope is small - yields a quantitative 
requirement on the tightness of £ and TZ. 

4. Singular Hopf normal form 

In a slow-fast system, an equilibrium point may cross a fold of the critical manifold. If it 
undergoes a Hopf bifurcation at 0(e) distance from the fold both in parameter and phase space, 
we follow [H] and refer to this as a singular Hopf bifurcation. Singular Hopf bifurcation occurs in 
generic one parameter families of slow-fast systems. 

We use a normal form for singular Hopf bifurcation in systems with one fast and two slow vari- 
ables proposed by Guckenheimer [9j as an example system for the computations of slow manifolds 
presented in this paper. The normal form is given by 

ex = (y — x 2 ) 

(8) y = z — x 

z = —fi — ax — by — cz, 

which depends upon the four parameters /i, a, b, c as well as e. An er-dependent scaling transfor- 
mation eliminates e as a parameter: set 

(9) (X, Y, Z, T) = (e~ 1/2 x, e~ l y, e~ 1/2 z, e~ 1/2 t) and (A, B, C) = (e 1/2 a, eb, e 1/2 c) 
to obtain 

X' = Y - X 2 

(10) Y' = Z-X 

Z' = —fx — AX — BY — CZ 

Guckenheimer [H] studies invariant manifolds in the phase space of system ffTD]) at different 
system parameters: The branch of the critical manifold y = x 2 where x < perturbs into a 
repelling slow manifold, while the branch where x > perturbs into an attracting slow manifold. 
In large regions of parameter space, an equilibrium that has undergone singular Hopf bifurcation 
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FIGURE 5. Phase space of system (TIP]) when (fx,A,B,C) = 
(0.0015709,-0.05,0.001,0.1). A repelling sheet of the slow manifold is plot- 
ted in dark blue, an attracting sheet in cyan. A selection of trajectories in the 
unstable manifold of the equilibrium near the origin are plotted in red, a strip of 
unstable manifold that escapes from the fold region is shaded in magenta. 

is a saddle-focus with a two-dimensional unstable manifold that is initially bounded by a periodic 
orbit. As the parameter fi is varied, the unstable manifold grows and eventually intersects the 
repelling slow manifold, first tangentially and then transversally. In the following, we refer to such 
a tangential intersection of the equilibrium's unstable manifold with the repelling slow manifold 
as a tangency or tangency of invariant manifolds. Figure [5] shows selected objects in the phase 
space of system ( llOp at a parameter just after the tangency. The tangency is a codimension 1 
phenomenon. Note that since slow manifolds are not unique, there is some ambiguity in what it 
means for a slow manifold to intersect another manifold. We will introduce a definition to deal 
with this in section [7J 

The tangency bifurcation is evident in the organization of phase space by invariant manifolds, as 
it separates regions in parameter space where trajectories in the unstable manifold of the singular 
Hopf equilibrium have different possible limit sets: after the tangency, trajectories can escape from 
the fold region, whereas, before the tangency, the unstable manifold is confined to the fold region. 
This change is significant in many other slow-fast systems: suppose that a system with one fast 
and two slow variables has an S-shaped critical manifold as well as a singular Hopf bifurcation 
followed by a tangency. The system introduced by Koper [20l [22| [23| [5] is an example. Figure [6] 
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(left pane) shows a trajectory in Koper's system just after the tangency bifurcation, together with 
the position of the critical manifold: the trajectory starts in the vicinity of the singular Hopf 
equilibrium and goes through a spiraling motion until it leaves the fold region to the left. It lands 
close to the critical manifold on an attracting slow manifold, which it follows to a fold before 
"jumping" to another attracting slow manifold and returning along this slow manifold back to the 
vicinity of the singular Hopf equilibrium point. The process now repeats, leading to a time series 
like the one shown on the right pane of Figure |6j Such patterns with alternating large-amplitude 
and small- amplitude oscillations are known as mixed mode oscillations in the literature [5]. 




x t 

FIGURE 6. Mixed mode oscillations in Koper's system with system parameters e± = 
0.1, 62 = l,k = —10, A = —7.50: the left panel shows a trajectory in the xy-plane, 
where x is a fast variable while y is slow. The singular Hopf equilibrium is marked 
with a black dot, the critical manifold is drawn with a dashed black line. The right 
panel shows the time-series of the x-coordinate of the same trajectory. 

Guckenheimer and Meerkamp [TT] present a detailed analysis of the bifurcation structure of 
the singular Hopf normal form. This includes extensive numerical results on the position of the 
tangency bifurcation in the five-dimensional parameter space of the normal form. The position 
of the tangency curve was computed in two ways: (1) using the numerical continuation software 
AUTO [39] and (2) using custom MATLAB code. In method (1), a boundary value problem is 
set up to track a trajectory segment that starts on a fixed ray in the linear approximation of the 
unstable manifold of the singular Hopf equilibrium point and follows the repelling slow manifold 
for a substantial period of time. The latter is achieved by requiring it to have a sufficiently long 
time-length and ending on the parabola y = x 2 + 5, thus forcing the trajectory to remain very 
close to the repelling slow manifold nearly to the end of the trajectory segment. A tangency of 
invariant manifolds corresponds to a fold of a solution to this boundary value problem, i.e., a point 
where two solutions approach each other and vanish together as the parameter is varied. Such 
folds can be detected by AUTO in a one-parameter continuation and continued in two parameters, 
thus enabling detection and continuation of the tangency bifurcation. Method (2) is based on the 
heuristic that the tangency separates regions in parameter space where trajectories in the unstable 
manifold escape the fold region from regions where trajectories in the unstable manifold remain in 
the fold region. A grid of initial conditions in an approximation of a fundamental domain of the 
unstable manifold is integrated numerically for a sufficiently long time. If the sample trajectories 
in the unstable manifold limit to an attracting periodic orbit or another bounded attractor, the 
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parameter is before the tangency. If at least one trajectory leaves the fold region, the parameter 
is after the tangency. Repeated applications of the above steps in an interval bisection method 
determine an approximate position of the tangency in parameter space. Note that neither method 
(1) nor (2) is rigorous. In particular, neither of the two methods establishes bounds for the position 
of the repelling slow manifold or of the unstable manifold of the equilibrium. Section [7] presents a 
rigorous method to compute a tangency. 



5. Detailed description of the method for the singular Hopf normal form 

In this section we give a detailed description of our method for computing enclosures of slow 
manifolds, applying it to the system from Section H] as an example. Most of the details generalize 
to any system of the form (j3j). In the description, we comment on nontrivial differences between 
the general case and the example at hand. 



5.1. Constructing the triangulation. The first step of our algorithm is to triangulate a portion 
of the critical manifold So- On a normally hyperbolic piece of the critical manifold, d x f ^ 0. The 
implicit function theorem implies there is locally a function ho(y,z), such that / (ho(y , z) , y , z) = 
0. In the singular Hopf normal form, h is given explicitly as h^[y,z) = ±y/y with domain 
D = [y m ,yM] x [Zmi z m\ C M 2 . For other systems, any suitable method for finding a sufficiently 
accurate approximation to ho{y,z) can be used. 

To construct the vertices of a Delaunay triangulation of So, as shown in Figure QJa), we start 
with a triangulation of the domain of ho, but want the diameter of the triangles on Sq to be almost 
uniform. Setting n{y,z) = \\Vh\\, k = y/ (y M - y m ) 2 + (z M ~ z m ) 2 /d, and 



k(y,z) :-- 



k 



1 + n(y,z) 



with d G Z + to be chosen later, we select the following points in the (y, z) plane as vertices of a 
triangulation: 



'ID 



(12) 



(Vo, z ) 
(Vi, zo) 
(yu zo) 



+ k{y i ^ 1 ,zo),Zo), 
(yM,z m ), 
{yi,Zj_ x + fc(2/i,Zj-_i)), 

(yi,z M ), 

< % < I, < < J, 



if < yi-i + zo) < y M , 

if < y M < y%-\ + k(yi-i, z ), 

if zj-! < Zj_i + k(y h Zj^i) < z M , 

if zj-x < z M < Zj-i + k(yi, Zj-t), 



Note that these points are aligned along lines parallel to the fold curve x = y 
Let T denote the Delaunay triangulation generated by the set 

{(Vi,*) : <i < l,0<j(i) < Ji}, 



where d x f = 0. 



and K-o its lift to So, using the map tTq 1 : (y, z) \— > (ho(y, z), y, z). Clearly n^ 1 is a homeomorphism; 
i.e., the set of vertices, edges, and faces of /Co, denoted by V(JCo), E(JCq), and F(JCo), are defined 
by tTq 1 (V (T)) , 7T _1 (£'(T)), and 7r " 1 (F(T)), respectively. T and /Co are shown in Figures [TJ^ a) and 
fjjb), respectively. 
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5.2. Constructing perturbed triangulations. Our next step is to perturb /C , as illustrated 
in Figure [21 so that it lies closer to the slow manifold S £ we are trying to enclose. Fenichel theory, 
|18| ; guarantees that for e > sufficiently small, S s is the graph of a function h e (y, z) with domain 
D and h £ (y,z) — ho(y,z) = 0(e). To compute triangulations, K. £ that approximate S e , we write 
h £ in the form 

he(y, z) = h (y, z) + eh^y, z). 
Substituting into the equation ex £ = f(h £ (y, z),y, z), we get that: 

f(h (y, z) + eh x (y, z) + 0(e 2 ), y, z)je = d y (h (y, z) + ehx(y, z))y + d z (h (y, z) + eh^y, z))z 

= d y h (y,z)y + d z h (y,z)z + 0(e) 

= d y h (y,z)g y (h (y,z),y,z) 

(13) +d z h (y, z)g z (h (y, z), y, z) + 0(e) 

To compute d y ho and d z ho, we use that f(ho(y, z), y, z) = 0, and hence 

a h (., 7 \ _ d y f(h (y,z),y,z) 

o y rio{y, z — — r r, 

d x f(h (y,z),y,z) 

and 

d z f(h (y,z),y,z) 



d z h (y, z) 



d x f(h Q (y,z),y,z)' 
In addition, since f(h (y,z),y,z) = 0, 

(14) f(h £ (y,z),y,z) = ed x f(h (y,z),y,z)hi(y,z) + O(e 2 ). 

Thus, we can solve equation (I14p for h±(y, z), up to 0(e), and substitute for f(h £ (y, z),y, z) using 
(fT3|) . obtaining 

, / \ d y f(h (y,z),y,z)g y (h (y,z),y,z) + d z f(h (y,z),y,z)g z (h (y,z),y,z) 

hi(y, z) = — 2 1" 0(e), 

(d x f(h (y,z),y,z)) 

which in our case, considering h + (y, z) reads: 

(15) ht(y,z) = ^^. 
For h~(y, z), that we will use in Section [TJ we get: 

(16) h^y,z) = ^_Zl. 

We put 7TT 1 : (y, z) \-> (h (y, z) + ehi(y, z), y, z), and define: 

IC £ := n £ l o 7T (/C ). 

K. £ is our approximation to the slow manifold, shown together with So in Figure [2^b). Heuristically, 
it is 0(e 2 ) to S £ at the vertex points. 

Let er c denote the following map that moves points parallel to the x-axis: 

(17) a c : (x,y,z) ^ (x + cmax ({h^y, z)\, ^) ,y,z). 
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We define our candidate enclosing surfaces as: 

(18) C £jN := o\_ e/A r(/C e ) 

(19) TZ £ , N := a £/N (IC £ ), 



where N G R+. The initial choice for N in our implementation was N = 64, but we would have 
chosen a smaller N if that had failed. The verification step of the algorithm includes a loop that 
divides iV by a factor 2 upon failure and repeats the transversality test. Note that the region 
that is enclosed by L £) n and TZ £i n is disjoint from the critical manifold so long as N > 1. The 
construction of Sq, S £ , C £> n and TZ £ ,n is shown in Figure [2j 

5.3. Verifying the enclosure property. To prove that a slow manifold is located between 
£ £i n and 1Z Ej n, it suffices to prove that the vector field (jSJ) is transversal to each face of the 
triangulations, with opposite crossing directions for C £t N and TZ £ ^. For the remainder of this 
subsection, we restrict our attention to a single triangle. Local transversality, i.e., the verification 
of transversality on each face in the triangulation implies global transversality of C £ ^n and TZ £ .n- 
Let T be one face in C £ ,n or TZ £j n- We denote its vertices by V\, v 2 , and v% and its edges by 
Ci2; Ci3, and e 23 with the edge between the vertices Vi and Vj. To verify that the vector field 
is transverse, it suffices to prove that the inner product between the normal of the face and the 
vector field is non-zero. Note that in contrast to most work on slow-fast systems, this condition, 
which is the main condition checked by our algorithm, becomes easier to verify as e — > 0. The 
reason is that as e — > 0, the condition becomes essentially one-dimensional. We denote the normal 
to the face, normalized so that the first component is positive, by n(T). This is possible because 
the first component is zero exactly at the folds, where the critical manifold fails to be normally 
hyperbolic. With this notation, the condition that we have to verify is 

(20) F(x, y, z) ■ n{T) ^ 0, for all (x, y, z) G T. 
Condition (120)) is equivalent to a verification that 

(21) F(X 1 v 1 + X 2 v 2 + X 3 v 3 )-n(T)^0 for aU A* G [0, 1], X 1 + A 2 + A 3 = 1, 

which is an enclosure of the range of a function on a compact domain. This problem is the one we 
solve with interval analysis. Directly enclosing (12~T1) using interval analysis in order to verify that 
the function is non-zero is, however, not optimal. The reason is that the problem is sufficiently 
sensitive that we would have to split the A» domains into a very fine subdivision, and since this 
has to be done on each face, such a procedure would be prohibitively slow. 

Our actual approach is based on monotonicity; first we prove that F ■ n is monotone on the face 
and on its restriction to the edges. Then we compute F(vi) ■ n for the three vertices and verify 
that the interval hull of the results, i.e., the smallest representable interval containing the results, 
does not contain 0. Note that this amounts to showing that the dot-product does not change sign 
on the face. We introduce 

(22) G:=V(F-n). 

If G 7^ (0, 0, 0) on all of T then F ■ n has no critical points inside of T and we can restrict our 
attention to the edges, i.e. the boundary of T. Consider an edge e^- = {(1 — X)vi+Xvj : AG [0,1]}, 
and denote its parametrization by r(A). The scalar product F ■ n is monotone on the edge if 

0?£{F(r{\))-n) = G-{v J -v i ). 
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Hence, we arrive at the monotonicity requirements, which for the case at hand are much easier 
to verify than (|2lj) : 



(23) 


(0,0,0) s 


t G(T) 




(24) 


s 


t G(e 12 )- 


(v 2 - v x 


(25) 


S 


t G(e 13 )- 




(26) 


S 


t G(e23). 


• (v 3 - v 2 



If the conditions ( 1231126"]) are satisfied we compute 

(27) F(vi) ■ n U F(v 2 ) ■ n U F(v 3 ) ■ n, 

where U denotes the interval hull. If (|27p does not contain zero, then the vector field is transversal 
to the face T. If f[2~3"j) holds but one or more of f l2"4"]|2"6"|) do not hold, then we add the appropriate 
F(eij) ■ n terms to (I27|) . 

5.4. Improving the bounds. If the previous steps of the algorithm are successful, they yield 
two surfaces C e ^ and TZ £ ,n, that have been proven to enclose the part of the slow manifold that 
is above [y m ,yM} X [z m , zm] m the (y,z) plane. Since N is fixed after the verification step we 
henceforth drop the indices on £ and 1Z. Our aim is to produce enclosures that are as tight as 
possible, given the mesh size. We, therefore, try to improve the enclosure. The procedure is 
illustrated in Figure |H 

We do this by iteratively updating each of the vertices in the triangulation by moving them 
towards each other along the segment joining them. This segment is parallel to the x-axis due 
to our earlier constructions. The moves are done in two steps: (1) a tentative move is made of 
a vertex, and (2) the transversality conditions of all faces attached to this vertex are verified. 
When the transversality holds, the vertex is fixed at its new position and we proceed to the next 
vertex. The efficiency of this procedure will depend on several factors, primarily the ordering of 
the vertices and how much the vertices are moved. By moving a vertex only a fraction of what 
seems to be possible, the effect of the ordering of the vertices can be minimized. The penalty of 
smaller updates is that the procedure has to be run more times. Larger moves might be possible 
if an appropriate sorting algorithm were used, but we have not found an effective and efficient 
sorting criterion. Instead, we heuristically determine an update factor that optimizes the accuracy 
vs complexity. Given a right vertex, Vr, and a left vertex, vl, such that tto(vr) = Ko(vl), we move 
each of them towards each other by an amount 

(28) ~ ll(^-^)||. 

We run the procedure to refine the enclosures of the slow-manifold several times, until no 
further improvement is possible. The quantity we use to measure the quality of the enclosures is 
the average distance between the two triangulations at the vertices. Let i denote the number of 
vertices of the triangulations; by construction C and TZ have the same number of vertices, edges, 
and faces. The only difference between £ and TZ is the values of the ^-coordinates. We put 

(29) rj(C,n) = - r \\v R -v L \\. 

v 1 

If the triangulation is fine enough 7] will be 0(e 2 ). This fact is investigated numerically in Section 
El 
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5.5. Cone fields. In order to ensure that there are manifolds inside of the set C enclosed by C 
and 1Z, we need to have invariant cone fields on C, as introduced in Section [3J In this subsection 
we describe how such cone fields - one horizontal and one vertical - are constructed. Recall, see 
[19] . that a standard horizontal or vertical cone for a phase space with variables (x,y) is a set 
{7 1 1 x\ I > ||?/ 1|} or {7 1| y || > ||x||}, respectively, and that a cone is the image of a standard cone 
under an invertible linear map. Equivalently, a cone is the set of points where a non-degenerate 
indefinite quadratic form is non-negative. Since horizontal and vertical cones are traditionally in 
the expanding and contracting directions, respectively, we will call the cone in the normal direction 
the vertical cone, and the cone in the direction of the slow manifold the horizontal cone. Also 
recall that a cone field is invariant if it is mapped into itself by the derivative of the dynamics, 
i.e., if the set where the quadratic form is non- negative is mapped by the derivative into the set 
where the quadratic form at the image point under the map is non-negative. 

For the case at hand we will use 7 = 1 for both the horizontal and vertical cones in an appropriate 
coordinate system, such that the normal direction is in the vertical cone. A cone field is a map that 
associates a cone to each point of its domain. Given that f fTUj) only has one nonlinear component, 
we will use constant cone fields. To prove that the cone fields are invariant, we solve the variational 
equation for the time 0.0004 flow map, and use the eigendirections of the derivative of the flow as 
a basis, in which we represent the standard horizontal and vertical cones with 7 = 1. We verify 
that the vertical and horizontal cone fields are invariant, and that the vertical cone contains the 
fast direction, which ensures that S e defined in ([7]) projects injectively onto the slow variables, 
and, thus, is a graph over them. The flow time needs to be large enough for us to be able to prove 
the separation of the horizontal and vertical directions, but small enough that we do not move 
away too far in phase space. The value 0.0004 turned out to be a good choice. 

5.6. Algorithms. An implementation [43] of the method described above has been made using 
the IntLab package |41| for interval arithmetic. A detailed description of the main algorithm is 
given as Algorithm [TJ The algorithm that checks if the vector field is transversal to a face is given 
as Algorithm [2j Algorithm [1] takes a triangulation as input. That triangulation can be computed 
with any method, not necessarily the one outlined in Section 15.11 In Algorithm [2] the function 
sign(x) returns if G x. 

6. Numerical Results 

In this section we describe the results of several experiments illustrating the behavior of the 
enclosure computations. Given a system and a domain, there are two numbers that can be 
changed, the number d, which controls the mesh size, and the value of e. In the experiments 
below, we use the normal form, ([S]), for the singular Hopf bifurcation discussed in Section 3. We 
choose the same values of the constants as in the first part of /i = 10~ 2 , A = —0.05, B = 0.001, 
and C = 0.1. We enclose the branch of the critical manifold {y = x 2 } with x > 0. The results of 
four experiments are described below, in each of them we present the results as a plot of 77 vs e. 
In the first experiment, we fix the domain as a small strip: y G [0.01,0.2], z G [—0.01,0.01] and 
give the results for several values of 1 (defined implicitly by changing d). In the second, we take 
a square domain: y G [0.01,0.2], z G [—0.095,0.095] for comparison. Our third example analyzes 
the effect and usefulness of the tightening step described in Section I6~3"l In our fourth example, we 
investigate the heuristic constant 8 in the denominator of (|2"8"|) : the domain and constants are from 
the first example with its finest mesh. Note that our domains are such that y < 0, which means 
that the assumptions from Section [3] are satisfied, i.e., all trajectories with initial conditions in C 
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Algorithm 1: Implementation of the main algorithm 



Data: (f,g y ,g z ), h ,T, e 
Result: L, 1Z, n 
1 forall the (y, z) € T do 

o U („. „\ _ dyf(h (y,z),y,z)gy(h (y,z),y,z)+dzf(h (y,z),y,z)g z (h (y,z),y,z) , 

3 end 

4 iV = 64; 

5 transversal=false; 

6 iVi* 1 = T -number Of Faces; 

7 while -^transversal & N > 2~ 18 do 

8 x ie/t = h (y,z) + hi(y,z) - e/N\hi(y, z)\; 

9 bright = h (y,z) + hi(y,z) + e/N\h 1 (y, z)\; 

10 if getTransversality(T, xi e f t ) = —getTransversality(T,x r i g ht) = NF then 

11 transversal=true; 

12 else 

13 N = N/2; 

14 end 

15 end 

16 if -^transversal then 

17 exit(FAIL); 

18 end 

19 r/ = 1; 

20 r) new = 0; 

21 while rj new < n do 

„ ll x ie/t~ x righil . 

•Kleft — •Klefti %right — bright j 

forall the 1 < i < t do 

irz = T.adjacentFaces(i); 
xieft(i) = xufS) + Q.l2h{x right {i) - x left (i)); 

if getTransversality(tri,xi e ft,T.y,T.z) = —getTransversality(tri,x r i g ht,T.y,T.z) 
tri.numberOf Faces then 

Xleft{i) = Xl e ft(i); 



22 
23 
24 
25 
26 
27 



28 
29 
30 
31 
32 
33 

34 
35 
36 
37 
38 



39 



else 

Xleft(i) = xi e ft(i); 
end 

x r ight(i) = x righ t(i) - 0.125(x right (i) - xi e ft(i)); 

if getTransversality(tri, xi e ft,T.y,T.z) = —getTransversality(tri,x r i g ht,T.y,T.z) 
tri.numberOf Faces then 



else 



x 



right \ 



X 



right<$)i 



end 



end 



W^left •^rightW , 

ijnew — y/T~i ' 

40 end 

41 C = Triangulate(T .Triangulation, xi e ft, T.y, T.z); 

42 1Z = Triangulate(T .Triangulation, x rig h t , T.y, T.z) 
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Algorithm 2: get Transversality(Triangulation, Vertices) 



Data: F = (f,g y ,g z ), T(Triangulation,Vertices) 

Result: Intersections 

NF = T .number Of Faces; 

Intersections = 0; 

forall the 1< i < NF do 



4 


n 


= T-Normal(i); 


r\ 
tj 


{vi,V2,v s ) = T.Vertices(i); 


6 


(ei2, 


ei3, 


e23) = T.Edges(i); 


rr 
I 


G = 


V(F(T.Face(z)) • n); 


a 
o 


if G G then 


9 




Intersections^ = sign(F(T .Face(i)) ■ n); 


10 


else 






11 




G\2 - 


= V(F(e 12 ) • n) • ei 2 , G 13 = V(F(e 13 ) • 


12 




if $ 


t G12G13G23 then 


13 






Inter section s+ = sign{F{v\) ■ n U ^(^2 


14 




else 




15 






forall the a £ {12, 13, 23} do 


16 








if E G a then 


17 








| F a = F(e a )-n 


18 








else 


19 








| F a = F{v ai )-nUF{v a2 )-n 


20 








end 


21 






end 


22 






Intersections^ = sign{F\2 U F13 U F23) 


23 




end 




24 


end 






25 end 









?23 = V(F(e 23 ) • n) ■ e 23 ; 



leave in both forward and backward time, and tangencies of the vector field with DC occur along 
a plane where they have quadratic tangency. 

During the computations we use the function G defined in (1221) to prove the monotonicity 
properties that enables us to efficiently prove transversality. We note that for the example at 
hand, G is 

- 2 fn x -n y + ^n z 
_ oi 

V y/e Uz 

A trivial calculation shows that G = (0, 0, 0) if and only if x = — 25-y/e and n is a multiple of 
(1, ^7=, 1000), so monotonicity always holds on the right branch of the critical manifold. 

6.1. Varying 1. The convergence rate of the enclosures at the vertex points should ideally be 
0(e 2 ), since we have corrected for the linear term in the asymptotic expansion of h e . Our interpo- 
lating surfaces between the vertex points are, however, linear. The discretization size thus puts a 
curvature dependent restriction on the tightness of the enclosure. In Figure [TJ^a), we illustrate how 
77, for different values of 1 first decreases, but then reaches a plateau. Looking at r/ as a function 
of e, we see that as the mesh size decreases (a increases), 77 is approximately proportional to e 2 , 
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as expected. This gives a heuristic picture of how r\ depends on e: first, there will be a period of 
quadratic convergence, where the accuracy depends on e; while at the end, the accuracy oscillates 
around some fixed value and depends on the mesh size. In the intermediate region, the accuracy 
depends both on the ratio of time scales and the mesh size. In this region, the exponent will 
decrease from 2 to 0. Figure [7(b) illustrates the quadratic convergence region for the finest mesh 
size from Figure [7(a). 

As the plateau is reached s(e) defined in starts to increase. For e — 0.1 the enclosure is 
too wide for all trajectories inside to be slow. In Table [1] we give the slopes on the e interval 
[1CT 1 , 10~ 4 ] and bounds on the intervals where *J\e)s(e) < 1, for the various i values from Figure 
[7(a). We are only able to prove that the cone fields are invariant for e < 10 _L94 , which means 
that for e > 10 -1 ' 94 the normal hyperbolicity is too weak for the algorithm to work. Thus, for the 
finest mesh size, we prove that the computable slow manifold exists for 1CT 6 < e < 1CT 1 ' 94 . Finer 
meshes would prove the existence for smaller values of e. 




FIGURE 7. (a) log 10 rj vs — log 10 e for the various values of t specified in Table 
[TJ (b) Zoom in on — log 10 e G [1,4] for the value l = 162190. The least squares 
approximation of the slope in the steepest part (— log 10 £ G [2,3.5]) is —2.14, on the 
whole interval [1,4] it is —1.89. 



1 


1200 


4662 


18236 


40805 


72239 


112736 


162190 


Slope 


-1.40 


-1.58 


-1.70 


-1.76 


-1.82 


-1.86 


-1.89 


max — log 10 e 


4 


4.5 


5 


5 


6 


6 


6 



TABLE 1. The second row is the least squares approximations of the slopes of 
log 10 r/(— log 10 e) on the domain — log 10 £ G [1,4], for some different values of t. 
The third row gives the maximum value of — log 10 £, where the flow is slow, i.e., 



6.2. Larger domain. In this subsection, we redo the experiment above for a square domain. 
There are roughly the same number of triangles in the y and z directions, rather than having only 
a couple of faces in each {y = const} slice as we had in Subsection 16.11 The resulting n vs e graph 
is given as Figure [HJ We see that the results correspond to the coarser meshes in Figure [7(a), which 
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is natural, since a larger domain would require a larger number of faces. This illustrates that the 
results in Subsection 16. II do not depend on the specific thin slice in the z-direction that we chose to 
study. For the two discretization sizes in Figure El we have y/{e)s(e) < 1 for e > 10~ 4 and e > 1(T 5 , 
respectively. We are only able to prove that the cone fields are invariant for e < 1CT 2 ' 09 , which 
means that for e > 10~ 2 ' 09 the normal hyperbolicity is too weak for the algorithm to work. Thus, 
for the finest mesh size, we prove that the computable slow manifold exists for 10 -5 < e < 10 -2 09 . 




Figure 8. log 10 ?7 vs -log 10 e for the values i = 21810 and i = 194396. The least 
squares approximations of the slopes on the interval [2,4] are —1.27 and —1.52, 
respectively. 



6.3. The effect of the tightening step. The tightening step is the slow part of the algorithm, 
and our program spends the vast majority of its computing time performing this step. It is therefore 
interesting to see how the results of a fast version of the algorithm, without the tightening step 
compares, performance wise. We run the example from Section 16. 1\ with the highest precision 
(l = 162190), and compare the results. The 77 vs e graph of the results is given as Figure |HJ In 
this example, the program spends 92.7% of the computing time performing the tightening step. 
The total computing time in this case was 1526 seconds on a 3.2 GHz Dual-Core AMD Opteron. 
For the example at hand, it might not be worth the extra effort to compute the tightening step 
or all applications. We do need it, however, for the application in Section [71 

6.4. Varying the improvement rate. Our method contains a choice of the heuristic constant 
in the denominator of equation (128]) that regulates the aggressiveness of the tightening step. In 
this subsection, we present a study on how the results depend on this choice. We use the same 
model as above, the domain from Subsection 16. 1[ and the finest mesh size from Subsection 16.11 
- 1 = 162190. For the purpose of this study, we denote the denominator of equation fT28|) . by /. 
In Figure [TUl we display the results for / = 4,6,8. For larger values of I, the results are virtually 
indistinguishable from the I = 8 case. Typically, the updates mostly occur for smaller values 
of e. The reason is that for sufficiently small values of e the vector field is almost equal to the 
layer equation, which makes the transversality condition almost trivial. Therefore, less smooth 
triangulations will still work, and the updates will not violate the transversality conditions. 
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- log 10 e 

Figure 9. log 10 ?7 vs — log 10 £ for the value i = 162190, with and without the 
tightening step of the algorithm. 




- l°gio e 

Figure 10. log 10 77 vs — log 10 e for the value 1 = 162190, for updates with || {vr—vl) || 
divided by 4, 6, and 8. 

7. Tangencies 

In this section we give a proof that the singular Hopf normal form, given by fflQ|) . used here 
with (A, B, C) = (—0.07,0.001,0.16), undergoes a tangency bifurcation of the unstable manifold 
of the saddle equilibrium, and the repelling slow manifold. We will often refer to these manifolds 
as the unstable manifold and the slow manifold denoted by and ST, respectively. With a slow 
manifold for the rescaled system, we mean the image of a computable slow manifold for some e 
under the map (jUJ). 

Recall that (computable) slow manifolds are not unique. We therefore need to define what we 
mean by tangency, since if one choice of computable slow manifold is tangential, there will be 
other choices where the intersection is transversal. The natural setting is therefore to define when 
a one parameter family of slow manifolds is tangential to another manifold or family of manifolds. 

Definition 7.1. A smooth one parameter family of manifolds, {M^}^^^^, intersects a one 
parameter family of families of computable slow manifolds {C^j^e^o.Mi] tangentially if for each 
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choice of a smooth one parameter family of computable slow manifolds {S , At } (Ue [ At0)/11 ], G C M , 
there is a value of \i G (yU , /ii) snc/i i/iai S M and M M intersect tangentially. 

In our proof, we compute one enclosing region C that satisfies the requirements from Section 
[3] for all values of the parameter fi that appear in the proof. However, the computable slow 
manifolds might change with the parameter, since they are defined using ([7]). We prove that the 
one parameter family of unstable manifolds W^p^) moves through this fixed enclosing region, 
and that, as the family passes through C, it always has to have a tangential intersection with at 
least one of the computable slow manifolds, regardless of how the smooth one parameter family 
of computable slow manifolds inside of C was chosen. 

Theorem 7.2. For < e < 10~ 3 , the singular Hopf normal form undergoes a tangential 
bifurcation of a computable slow manifold and the unstable manifold of the equilibrium. The 
bifurcation occurs in the interval [fio,Hi] = [0.00454,0.004553] with fixed parameters (A,B,C) = 
(-0.07,0.001,0.16). 

The main argument in the proof of Theorem 17.21 is illustrated in Figure [TTJ We consider the 
intersections of and with a half-plane E. At /io, the two manifolds do not intersect each 
other in E. Notice that the unstable manifold seems to translate to the left relative to the repelling 
slow manifold as /i increases. At fii, the two manifolds intersect transversally in E. In the proof 
of the theorem, we formalize and prove these observations, and moreover show that the first 
intersection of the two manifolds is tangential. The vector field is transverse to E, so a tangential 
intersection of the manifolds in E corresponds to a tangential intersection in the 3-dimensional 
phase space. 

In the proof of Theorem 17.2} we will at times work with the singular Hopf normal form ([S]), 
and at other times with the rescaled singular Hopf normal form ([10]) . Recall that, in the rescaled 
system we use upper case variables and parameters (/i is scale independent). Note that we do not 
assert that the tangency of the manifolds is unique. We will first prove Theorem 17.21 for e = 10~ 3 . 
For smaller e, the result follows from the rescaling (Q and the following property of the relative 
slope condition for the singular Hopf normal form (recall that the tangency will occur in different 
parts of phase space for different values of e): 



2\X'\ 






+ 


Z'\ 



\fi\Y'\ + \ z ' 

since e > 0, \Y'\ > 0, and \Z'\ > 0, this is a non- decreasing function of e. 
Hence, if 

s(£) £ h" 

then 



\=Ve i s(e') < -^=Jls(e) < — =, for all < e < e. 



S[E 

y/e' \je \js' 

The existence of computable slow manifold at a particular value of e thus implies the existence of 
computable slow manifolds at all smaller values of e. Note that these computable slow manifolds 
will appear at different positions in the phase space for different values of e. 
Set-up. We will work with 

S := {(X, Y, Z) G M 3 : Z > -0.1693 + 0.16 (X + 1.353), Y = 2} 

and 

R := {(X, Y, Z) G E : -1.62 < X < -1.49, -0.169 < Z < -0.162}. 
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FIGURE 1 1 . An example family (thick line) and images of fundamental domains 
(solid curves) of for a selection of \i in [/i ,/ii], shown here intersected with E. 
The boundary of E is drawn as a dashed line, the rectangle R is drawn as a shaded 
region. As /i is varied in [/i , fii], the slow manifold only moves by amounts too small 
to be noticeable at the scale of the diagrams. For each included in the figure, 
we plot the first intersection of the trajectories with E, if the trajectory reaches E. 



We next list verifiable conditions that via Lemma 17.41 below will prove Theorem 17.21 Many of 
these conditions are illustrated in Figure [121 Let 

Ymi n ,Y max \ [fiQ, fli] —> R 

be continuous with Y min (fi) < Y max (n). Further define a 2-dimensional "box" by 

B := {(n,X,Y,Z) G [/i ,/ii] x : X = 7r x (p^,Y mi M < Y < Y max (fi)}. 

Note that the requirement (X,Y,Z) G uniquely defines Z as a function of (//, X, Y). Denote 
the corners of B corresponding to 

(fi, Y) G {(/xi, Y max (fii)) , (/i , Y max (fi )) , (fi , Y min (fi )), [fj,i, Y m i n {fii))} 

by {Mi, M 2 , M 3 , M 4 }. Denote the flow map of system (fTUj) from B to E, wherever it is defined, by 
ty. The next step of our construction is to introduce a number of assumptions, that are verifiable 
using validated numerics, i.e., they can be restated as a finite number of computable conditions. 
The geometry of these assumptions is illustrated in Figure [12J In Lemma [7.41 below we show that 
these assumptions are sufficient to prove Theorem 17.21 

Assumption 7.3. Assume that the following conditions are satisfied: 

(I) For fi G \po,/jLi], a family of repelling slow manifolds ST intersects R in a single family of 
curves C M that enters R at the top and exits R at the bottom. 



21 



JOHN GUCKENHEIMER, TOMAS JOHNSON, & PHILIPP MEERKAMP 



(II) The map ^ is defined on the three sides of B corresponding to Y = Y min (/j,), Y = Y max (fi) 
and fi = fio, and their images under \1/ lie in R and strictly to the right of ST H R. 

(III) The map ^ is defined on {(h, X,Y, Z) G B : h = Hi} and its image lies in R. Furthermore, 
ty(Mi) and \l/(Af 4 ) are strictly to the right of ST D R, and there exists a point M 5 G 
{(/x, X, Y, Z) G Bq : n = Hi} such that \P(M 5 ) lies strictly to the left of S^H R in E. 

(IV) The map ^ is well-defined on B Q . 



x 10 

-8.4 r -0.1635 




-9.8 1 1 1 ' -0.168 1 1 1 1 ' 1 

4.54 4.545 4.55 4.555 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 

(a) " *^ (b) 

FIGURE 12. Illustration of the assumptions made in Assumption 17.31 The box B Q 
shown in pane (a) maps into R C E as shown in pane (b). As fi is varied in [/x , Hi], 
the slow manifold (thick solid line) only moves by amounts too small to be noticeable 
at the scale of the diagrams. 

Lemma 7.4. Suppose that Assumptions 1 7. $ are satisfied. Then and intersect tangentially 
for some jj* G [/io, Hi] ■ 

Proof of Lemma \ 7.4\ Fix a family of slow repelling manifolds ST, \i G [fio, Hi]- Since all of B 
reaches E by Assumption 17. 31 IV. the existence and uniqueness theorem for ODEs implies that the 
map from B to E x [fio, Hi] is continuous. We may thus define the continuous function 

dist(/i, Z) := min {tt x (W;\ z D R) - 7T X (S;\ Z H R)) 

where Z is required to lie in the range of Z values of R and \z denotes restriction to Z. Consider 

H* = min{/i G [ho, Hi] '■ mindist(/i, Z) = 0}, 

z 

the existence of which follows from Assumptions 17.31 11 and l7.3I III. and the continuity of dist(/i, Z). 
Clearly ST* PI R and D R intersect in at least one point (X Q , Y Q , Z ). Moreover, 

itx(W^\ z OR)- 7T X (S M . \ z n R) > 

for the range of Z values that lie in R. Since and S^* are smooth surfaces in IR 3 , transverse 
to R, we can now consider the Taylor series expansion of 

Kx{w;*] z nR)-it x {S^\ z nR) 

at Xq with respect to X and conclude that its linear term must be zero. We have thus shown that 
the manifolds W™* and S M * intersect tangentially in R. □ 
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Subsection 17.11 below gives details on the verification of Assumption 17.31 1. Subsection 17.21 de- 
scribes in detail how E, Y min , and Y max are chosen, and provides details on the verification of 
Assumptions E3II, E3JIII, and QIV. 

7.1. Slow manifold computations. Showing that for fx E [fx , fxi], a family of repelling slow 
manifolds intersects R in a single family of curves is a straight-forward application of the 
methods developed in the earlier parts of this paper: we compute slow manifold enclosures for the 
rescaled singular Hopf system ( ITUl) over a domain that corresponds to 

1 < Y < 500, -0.169 <Z< -0.162, 

and for the singular Hopf parameter values 

{{fx, A, B, C) E M 4 : fx E [4.54, 4.553] x 10~ 3 , A = -0.07, B = 0.001, C = 0.16}. 

The actual computations for the enclosures are performed in the original singular Hopf coordinates 
of (jSJ), as described in earlier sections of this paper. Let eo = 10~ 3 , in the original coordinates of 
the singular Hopf normal form, (jSJ), the domain now corresponds to 



D \}Jmini Umax] ^ 



mini ^max] i 



where y min = 1.0 e , y ma x = 500.0 e , z min = -0.169 y/e^, z max = -0.162 v /e^, and the set of singular 
Hopf system parameters is 

[e, fx, a,b,c) EM 5 : e = e , fx E [4.54, 4.553] x 10~ 3 , a = -^1, b = °' 001 , c 



The enclosures obtained show that points (X, Y, Z) in the repelling slow manifold over D must 
satisfy 

-1.5726 < X < -1.5539. 

Moreover, the methods of Section [3] of this paper were used to check that at any parameter in the 
above-described set, ST is a graph over a domain D C So, and that s(sq) < 1.027. Again, note 
that the computation is independent of the choice of Eq, since a different choice of e would imply 
that we should enclose a different part of the phase space. Since the above z m i n and z max were 
chosen large and small enough, respectively, to conclude that the enclosed repelling slow manifolds 
enter R at the top and leave R at the bottom, we have shown the existence of the sought family 
of slow manifolds. 

Remark. Even though the slow manifold intersected with the section S in our case resembles a 
fixed straight line, enclosing it with the precision required for the proof to work is a hard problem. 
To determine rigorously the location of a slow manifold is difficult even in the easiest non-trivial 
cases. The problem is amplified in our case since we need high accuracy in the rescaled system, 

where the errors are blown up by a factor O (^-^ 

7.2. Unstable manifold computations. We now describe how Y min and Y max are chosen for 
Assumptions 17.31 11 and I7.3I III to be satisfied. Recall that Figure [TT| was obtained by examining 
trajectories in entire fundamental domains of for fx E [fXo, fXi], and that some of these trajecto- 
ries did not reach E. We chose Y min first and then Y max in such a way that B is on the one hand 
small enough for the map to E to be well-defined and its image to be in R, and on the other hand 
large enough for the images of marked points Mi, . . . , M 5 and of the boundaries of B to map to 
the left or right of S^ as required by Assumptions 17.31 11 and l7.3I III. 
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7.2.1. Computing Y min (fi). Let L C E be the line given by 

L := {(X, Y, Z) = (t, 3, -0.1678 + 0.16 (t + 1.535)) : t G R}. 

This line lies well within E and is parallel to <9E. It is moreover transverse to the parts that 
reach E, for all \x G \po,fii]. The boundary value problem (BVP) for the flow of the rescaled 
singular Hopf normal form with the following boundary conditions and fi as the continuation 
parameter is thus well-defined: 

• trajectories have to start in the unstable eigenspace of p M 

• trajectories have to start at an X coordinate equal to that of p^ 

• trajectories have to end on L. 

Note that there are multiple solutions to this BVP, as each trajectory in the unstable manifold 
satisfies the initial boundary condition multiple times as it spirals away from the equilibrium point 
Pp, but we choose one by selecting a fundamental domain for its endpoint near p^. The equilibrium 
points Pn = x 2 „, Xfj) satisfy the equation x M = —45 + \/A5 2 — 1000/i. We use a trajectory that 
initially has a Y coordinate approximately 10~ 4 larger than that of p^, deferring a discussion of 
the suitability of this distance to a remark at the end of this subsection. Solving the BVP with a 
shooting method, we find that the Y coordinates of the solutions to the boundary value problem 
are close to linear in \x on the interval [/xo, /xi]- We thus define 

Y min {fi) = x\ - 9.37888799540 x 10~ 5 + 0.640307054861539(/i - /i ), 

to be the linear function in \i that approximates the Y coordinates of the solution endpoints to 
the BVP. 

7.2.2. Computing Y max (fi). After inspecting diagrams similar to Figure [T2"} we defined Y max (fi) in 
an ad-hoc manner to be the linear function for which the box B contains 25% of a fundamental 
domain of WJf, and 40% of a fundamental domain of : 

Ymax(n) =%l- 9.628167607168 x 10~ 5 + 0.54980571 1513847(ii - /U ). 

7.2.3. Computing unstable manifolds. The complexity of the singular Hopf normal form makes 
it unfeasible to compute unstable manifolds analytically. We therefore begin by describing a 
method to rigorously compute the location of W%. We will use the method developed in Section 
[5] together with covering relations with cone conditions |36| and validated numerical integration 
|25[ [28| [29| [30] to enclose and propagate the manifolds, respectively. In our implementation [43] 
we use the software VNODE-LP |42| to integrate the system ( JTUj) . The computations are done 
using order 11 Taylor expansions in VNODE-LP. 

Since our proof relies heavily on the concept of h-sets and the method of covering relations we 
provide an informal introduction here. For a complete formal description of these concepts and 
methods we refer the reader to [37] [36] . In [37] h-sets and covering relations are introduced, and 
in [36| the concept of an h-set with cones is introduced together with the appropriate modification 
to the definition of a covering relation. An h-set is a compact hyperbolic like set, in the sense 
that it has expanding and contracting directions, in an appropriate coordinate system. An h-set 
is a set together with the coordinates. A map together with two h-sets, h±,h2, is said to satisfy 
covering relations if h\ is mapped across h 2 under the map. Across in this setting means that the 
boundaries of hi transversal to the expanding directions are mapped outside h 2 and the image 
of hi does not intersect the boundaries of h 2 transversal to the contracting directions. Using the 
Brouwer degree one can show, see [37J, that a cycle of h-sets with covering relations must contain a 
periodic orbit. An h-set with cones is an h-set together with a quadratic form Q, that describes a 
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uniform cone field on the h-set. The map is said to satisfy covering relations with cone conditions, 
if the quadratic form is increasing along orbits. Given recurrence, this yields uniqueness of periodic 
orbits. One can also use the cone conditions, see |36| . to prove the existence of invariant manifolds 
and propagate them along orbits, which is how they are used in this section. The bounds on the 
location of the invariant manifolds given by covering relations with cone conditions are Lipschitz. 
In particular around a fixed point one gets a cone, which bounds the location of the invariant 
manifold. The Lipschitz constant depends on the ratio of the positive and negative eigenvalues of 

Q. 

We construct an h-set with cones centered at p^ as a cylinder of size 10 4 and 10 5 in the (X, Y) 
and Z directions, respectively, with a cone with Lipschitz constant 0.1 defined by the quadratic 
form 





l 








Q = 





1 













-100 



We verify that covering relations and cone conditions hold for the time 6.3 map. This proves that 
the unstable manifold exists within the h-set, and yields an enclosure of the unstable manifold as 
a Lipschitz graph with Lipschitz constant 0.1 over the disc: 

{(X,Y) :||(X-x„F-xJ)||<10- 4 }. 

To further contract the enclosure for a given value of (X,Y), we partition the line segment over 
(X, Y) in the cone, and integrate backwards for 100 time units or until the trajectory passes the 
cone. Subsegments that leave the cone in backwards time are removed, and we use the interval 
hull of the remaining subsegments as our new bound of a point in the unstable manifold. The 
covering relations with cone conditions prove that each remaining subsegment over (X, Y) contains 
a unique Z value such that (X, Y, Z) G W". 

Given an initial enclosure of a point in we propagate it forwards by integrating ffTUj) until it 
hits E using VNODE-LP. To integrate the top and bottom of B , i.e., the boundaries of B where 
fj, is not constant, and the interior of B , we consider a 4 dimensional phase space by appending 
fi — to (fTU|) . This procedure stabilizes the numerical behavior of the propagation of the unstable 
manifold. 

7.2.4. Verifying Assumptions \7.$, (II-IV). Using the method described in Section [7.2.31 one can 
now subdivide OBq into small subsets, compute an interval enclosure of each subset, and use 
validated numerical integration to show that Assumptions 17.31 11. 17.3I III. and I7.3I IV are satisfied. 
In practice, this requires some experimentation: if the subsets are too large, wrapping effects 
in the numerical integration will make the verification of Assumptions 17.31 11. 17.3I IIL and I7.3I IV 
impossible. On the other hand, the computing time for the entire verification of Assumption 17. 3I II 
is approximately proportional to the number of subsets to be integrated numerically. The bounds 
on \[>(<9I?o) and \I/(Mj), for z = 1,4, and 5, are given in Table [2j 

Remark. Note that since the position of as well as the map to E are computed using 
interval arithmetic, their computed positions have errors due to over estimation associated with 
them. These errors have to be taken into account when choosing the Y value at which to place 
the half-plane E, the interval boundaries fi and and the functions Y min and Y max . Generally, 
placing E at greater values of Y results in tighter bounds for the slow manifold, and the repelling 
nature of the slow manifold spreads trajectories that were initially close in the fundamental domain 
far apart, making it easier to verify Assumptions l7.3l (II-IV). We found the size 2 x 10 -4 of the 
h-sets constructed in Section 17.2.31 to be large enough to keep the validated numerical integration 
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tf(0B o (y min )) 


#(d5 (^»a*)) 


X 


i c4fj8 


i 6102 
-"-•5156 


i c4<;2 
-"-• J 236 


1 cii68 
l.J 107 


z 


-0-16J5 


-0.16i 


-0.16- 


-0.16S 





^(Mi) 




V(M 5 ) 


(b) X 


1 . 5229 


-1.535? 


-1-57S 




—0.164^ 


-0.167? 


-0.166? 



TABLE 2. (a) The image of 8Bq under ^ . (b) The image of the marked points on 
the dBo(ni) line under \F All images are in the interior of R. The computations in 
(a) and (b) prove Assumptions 17.31 11 and l7.3I III. respectively. 



to E short enough to not accumulate prohibitively large errors, while being small enough to be 
efficiently computable. 

Remark. To give further insight into what happens after the bifurcation we note that the fol- 
lowing set is forward invariant. For other values of the parameters, similar sets can be constructed. 
For k > B, 

X<—^-, X*>(l + k)Y, Y> 1 -±^, \X\>\Z\. 
A + C k 

We verify that the above conditions are satisfied, with k = 2, for the point M 5 . Thus, X — > — oo 
and Y — > oo for a part of the unstable manifold past the tangential bifurcation. 

8. Summary and Discussion 

Computation of the slow manifolds in a normal form for singular Hopf bifurcation served as 
a case study for this paper. A singular Hopf bifurcation in slow-fast systems with two slow and 
one fast variable occurs when an equilibrium point crosses between attracting and repelling slow 
manifolds. The dynamics associated with this crossing - a folded saddle-node type II in the singular 
limit - is complicated. The small amplitude oscillations emanating from the equilibrium point are 
part of mixed mode oscillations in some examples, notably the model originally studied by Koper. 
Subsidiary bifurcations occur, including tangency between the repelling slow manifold and the 
two dimensional unstable manifold of the equilibrium point. Tangency bifurcations form part of 
the boundary of the parameter space region in which mixed mode oscillations occur in the Koper 
model, making them essential to understanding global aspects of the dynamics in this and other 
systems. Since there are no analytic methods for locating the tangency bifurcations, this paper 
uses verified computing methods to prove the existence of tangency bifurcations between a slow 
manifold and an unstable manifold of an equilibrium point for the first time. 

Some of our ideas generalize to the case of slow manifolds of saddle type. To compute normally 
hyperbolic manifolds of saddle type, see e.g. j3], one usually first computes the manifold's stable 
and unstable manifolds, and then intersects them. To compute a saddle slow manifold in a 
three dimensional ambient space using our ideas, one could compute enclosures of the stable 
and unstable manifolds, as presented in this paper. The existence argument given in Section 
|3] can be modified to this setting, under appropriate assumptions on the dynamics on the slow 
manifold. Generalization to slow manifolds of saddle type in higher dimensional ambient spaces 
is substantially more challenging. 

We made several design decisions while constructing our algorithm for computing slow manifolds. 
This section discusses details of several and motivates our choices. 
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Our enclosures were constructed as pairs of enclosing transversal piecewise linear surfaces. 
There are several alternative approaches to how to construct and refine the vertices of the 
enclosing triangulated surfaces C and 1Z. For the examples in Sections [6] and [7] we used 
rectangular patches in the domain of the slow variables. Instead, one could construct the 
triangulations of the original domain in the slow variables by considering a dynamically 
defined region, constructed by flowing a set of initial conditions on the critical manifold 
with the slow flow, and use a discretization of those trajectories as the vertices of the 
triangulation. 

We considered other possibilities for moving vertices in Section I5.4[ namely, to move them 
along trajectories of the flow of (jSJ), or to move them along the normal of the triangulation. 
Both of these methods have serious disadvantages. When moving vertices along the flow 
of the system, we have to carefully check whether the vertices are moved past edges, 
thereby destroying the integrity of the triangulation. If the triangulation remains a graph 
over the (y, z) domain, it is possible to generate a new triangulation by a Delaunay-type 
algorithm, and lift it to the surface, but if two vertices flow to the same (y, z) coordinate 
this is no longer possible. Additionally, this method of moving vertices moves the two 
enclosing surfaces by different amounts, so that we obtain an enclosure of a smaller part of 
the slow manifold. A third drawback is that the triangulations might develop very acute 
triangles. Finally, numerical integration of a large number of vertices is slow compared to 
the approach that we use. Moving vertices along the normals combines the worst of both 
methods: we no longer control the triangulations, and we might introduce violations of the 
transversality conditions. 

The tightening procedure described in Section [5741 only updates one vertex at the time, i.e., 
we move one vertex a big step and if all the faces attached to it are still transversal to the 
flow, then we move it. An alternative would be to move not only the vertex itself, but at the 
same time all vertices attached to it by an edge. Such a procedure would work as follows: 
when it is one vertex' "turn", only update it by a fraction of its potential improvement, and 
simultaneously move the ones it attaches to, by a smaller amount. The smaller neighbour 
updates should be such that the expected value of the total update of each vertex stays 
the same as in Section 15.41 The benefit of such an approach is that the triangulation is not 
skewed as much in each step, so it should be easier to verify the transversality condition. 
In practice, however, the gain of this approach is negligible, compared to a slight increase 
of the denominator of f )28|) . There are also disadvantages of such an approach, primarily in 
its computational complexity. Each time an update is made, one has to not only locate all 
its neighbouring vertices and update them, but also locate all of their neighbouring faces 
and check the transversality condition on them. In the results presented in Section [HI we 
thus only update one vertex at the time. 

We construct invariant cone fields on C to prove that it contains normally hyperbolic 
locally invariant manifolds. We constructed these manifolds by flowing a "ribbon" around 
the inflowing boundaries of the enclosure. The property that our enclosures were aligned 
with the flow in the sense that for one of the slow variables the vector field is non-zero, 
was crucial for proving the existence of computable slow manifolds. In general one could 
also use the invariant cone fields to show that the graph transform is well defined, by 
adapting the method in [19] . To prove the convergence of such a scheme would require very 
careful estimates of the expansion and contraction rates, and the norms of the nonlinear 
components of the vector field. An alternative is to define an extension of the vector field 
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outside of C that has a slow manifold that is invariant rather than just locally invariant. 
Global invariance together with normal hyperbolicity would give a unique manifold for 
the extension using the technique from [4j. Given normal hyperbolicity, ensured by the 
existence of the cone field, either method would give the existence of a (non unique) C 1 
normally hyperbolic manifold, which is the graph over the slow variables. Either of these 
approaches, however, include many subtle details that need to be clarified for the case at 
hand. 

• If the mesh size of piecewise linear enclosing surfaces remains fixed as e decreases, then the 
curvature of the slow manifold becomes a a limiting factor in the tightness of enclosures. 
With smoother enclosing manifolds, tighter enclosures are likely to be possible. We did 
not attempt this because the transversality calculations for piecewise linear systems were 
particularly simple in the singular Hopf normal form we studied. 
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